
clc
clear all
close all

Parameters;
iterb=1;
load('Parameters.mat');
betar = 0.8;
betaa = betar;


gg = (1+1.85./100)^(sig-1)-1;
BGP_allocations;
difff=10;


 while abs(difff)>10^(-10)

    
 
  %%Use definition of T = \sum A to update growth rate of T sector
     
     for i=1:M
         for j=1:M
Delta_eps(i,j) = ((eps(i,j))/(eps(i,j)+gg))*(1+gg)*lam(j)*(P(j)*H_r(j)/((Yw)))^betar;
       end
     end
     
[a,b]=eig(Delta_eps);
[c,vv]=max(diag(b));


%T=T./T(M);


LHS2 = gg;
RHS2 = max(real(eig(Delta_eps)));
difff = ([(LHS2'-RHS2')]');
scl=0.5;
T=real(a(:,vv));
Tnew=T./T(M);
betar = 0.85*betar+0.15*(betar  - scl*difff); 
%betar=betar(iterb);
betaa=betar;
iterb=iterb+1;  

BGP_allocations;


 end

 
 %%Productivity (from ACR formulation)
 
 for i=1:M
    Y_L(i)=Q(i)*(alpha(i))^alpha(i)*(1-alpha(i))^(1-alpha(i))*((1/mbar)*Tnew(i)^0.25/pi_shareT(i,i)^0.25)^(alpha(i))*((1/mbar)*TNT(i)^0.25/pi_shareNT(i,i)^0.25).^(1-alpha(i));
   
 end

 
 save 'Parameters_calibration' eta Bbar sig gg mbar alpha distT betar TNT lamNT L distT distNT LT LNT tau rho lam epsbar Q ipr chi ipr_rho gam_barg Ha w gg Tnew beta